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ABSTRACT 


A multi-dimensional numerical model has been' developed for 
the unsteady state osci] latory combustion of solid propellants 
subject to acoustic pressure disturbances. Including the gas 
phase unsteady effects, the assumption of uniform pressure across 
the flame zone, which has been conventionally used, is relaxed 
such that a higher frequency response in the long flame of a 
double-base propellant can be calculated. The formulation is 
based on a premixed, laminar flame with a one-step overall 
chemical reaction and the Arrhenius law of decomposition with no 

t 

condensed phase reaction. In a given geometry, the Galerkin 
finite element solution shows the strong resonance and damping 
effect at the lower frequencies, similar to the result of Denison 
and Baum. Extended studies deal with the higher frequency region 
where the pressure varies in the flame thickness. The nonlinear 
system behavior is investigated by carrying out the second order 
expansion in wave amplitude when the acoustic pressure 
oscillations are finite in amplitude. Offset in the burning rate 
shows a negative sign in the whole frequency region considered, 
and it verifies the experimental results of Price. Finally, the 
velocity coupling in the two-dimensional model is discussed. 


l 
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NOMENCLATURE 

a speed of sound 

B frequency factor for gas phase reaction 
C p specific heat at constant pressure for gas 
C s specific heat of solid 

E gas phase activation energy 

E s surface activation energy 

F eigenvector 

6 source term in numerical formulation 

h heat of combustion per unit mass 

H heat of reaction 

I order' of perturbation 

k thermal conductivity 

l characteristic flame length 

L latent heat of vaporization 

m mass flux 

M b Mach number 

n order of chemical reaction 

P pressure 

Pr Prandtl number 

r burning rate 

R gas constant 

R c dimensionless distance from surface in Eq. (22) 
t time 

T temperature 

u axial gas velocity 

Uj velocity vector 



u c core velocity in Eq. (22) 

v normal gas velocity 

w reaction rate 

x axial distance parallel to the surface 
X solution vector 

y normal distance from the surface 

Y mass fraction of fuel species 

z oxidizer-fuel ratio 

a thermal diffusivity 

6 ratio of solid to gas density 

s temperature exponent in Eq. (7) 

y specific heat ratio 

e perturbation parameter 

f = k 8 *C p */k*C s * 

\ eigenvalue 

\ y Lagrange multiplier 

t = k*/k s * 

p density 

u frequency 

Subscripts and Superscripts 
* dimensional quantity 

steady-state mean value 
a spatial variable 

(i) i th order perturbation coefficient 

+ gas side of interface 

- solid side of interface 

c propellant cold side 
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D time-dependent variable 

i,j vector quantity 

I time- independent variable 

0 mean value at flame edge 

s solid phase 

a, B finite element global node number 
t complex conjugate 


1 . INTRODUCTION 

An accurate analysis of combustion instability resulting 
from the oscillatory burning of solid propellants has long been 
of major concern. Such an analysis is characterized by the ^ 
acoustic admittance or response function, a proper measure of 
instability. However, experimental measurements of the 
coefficients are difficult; thus, it is desirable that analytical 
or numerical calculations be implemented whenever possible. 

Most of the past investigations on combustion instability 
have been centered around a one-dimensional quasi-steady analysis 
limited to pressure coupling [1-19], The representative studies 
were accomplished by Hart and McClure [2], Denison and Baum [4], 
and Culick [12-14]. Culick's review article [12] on homogeneous 
propellants summarizes the state-of-the-art up to the late 1960's 
and discusses the limitations of the analyses. Recently, several 
works concerning unsteady state problems have been attempted [20- 
27]. T'ien [21] introduces improved gas dynamics in order to 
recover the quasi-steady limitations by considering both gas and 
solid phases in an unsteady manner. However, the pressure is 
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still assumed to be constant and a function of time only, which 
has generally been used. Consequently, as indicated, in the case 
of the long flame of a double-base propellant, application of the 
model is limited to lower frequencies. 

Flandro [23,24] presents the first attempt to examine the 
response functions under the effect of incident acoustic waves in 
a two-dimensional analysis. A detailed formulation is extended 
to the second order perturbation system in terms of the Prandtl 
number and the Mach number. The effects of viscosity and heat 
transfer in the gas phase are included. The numerical results 
for the first order system appear to be comparable to those of 
T'ien [21]; however, the results of the second order response 
function are not clear. Moreover, details of matching conditions . 
are not given. The requirements for a more realistic model of 
combustion instability have lead to quite a number of rigorous 
studies. One recent work deals with the heterogeneous 
propellants summarized in Cohen's review paper [30-33]. 

The major discussions of the acoustic admittance or response 
function in the literature are concerned with pressure coupling 
usually applicable in the linear stability. Computations of the 
response function for velocity coupling are important in a 
nonlinear process; however, they remain in a state of infancy, 
although some initial attempts toward this subject have been made 
[34-38], Observations indicate that combustion oscillations are 
time-dependent and often nonlinear, as influenced by turbulent 
flow fields which may lead to erosive burning and unstable 
oscillations. Since the complicated physical phenomena cannot be 
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described analytically, some approximations and simplifications 
are still introduced to the theoretical formulations in the 
numerical analyses. Recently, Chung and Kim [26] introduced the 
effect of radiative heat transfer on the combustion instability 
in solid rocket motors based on the research of Chung and Kim 
[27], Further research has been reported concerning the 
calculation of response functions in multi-dimensional combustion 
phenomena [28]. The finite element method is introduced and 
complicated boundary conditions are handled easily by means of 
Lagrange multipliers [29]. 

On the other hand, much experimental effort has been devoted 
to understanding the oscillatory combustion mechanism using the f 
T-burner and L* -burner [39-42]. Extensive experiments may be 
found in the recent series of papers conducted by Levine and Baum 
[43,44]. Also, in the work of Caveny et al. [45], the 
oscillatory velocities in the solid propellant flames subject to 
pressure coupling are measured directly using the laser Doppler 
velocimetry instrumentation. 

The purpose of the present study is to examine the 
combustion instability induced by acoustic disturbances in the 
multi-dimensionally unsteady state in such a way that the upper 
limit of the acoustic wave frequency, above which an analysis 
with the assumption of the uniform pressure field in the flame 
zone cannot be applicable, is relaxed. Thus, very high pressure 
burnings (where any dynamic effects are to be important) and the 
long gas flame (such as burning of double-base solid propellants) 
could be possibly sought. 
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For small amplitude oscillations, the nonsteady governing 
equations are linearized by means of the first and second order 
perturbation expansions and solved numerically using the finite 
element method. The theoretical model is still based on a 
homogeneous propellant, Lewis number of unity, a second-order, 
single-step forward chemical reaction, and vaporization in an 
Arrhenius fashion, with no erosive burning. 


2 . ANALYSIS 

For the simplicity of the combustion modeling of solid 
propellants, the gaseous flame is assumed to be multi- 
dimensional, premixed laminar, and calorically perfect, and a ^ 
one-step forward chemical reaction occurs. The combustion of a 
solid propellant is approximated by the Arrhenius Law. The Lewis 
number is given as unity for the explicit relation between 
concentration and temperature. Thus, the conservation laws for 
the multi-component reactive system in the gas phase are 
represented as follows: 

Continuity 

dp 

— + (p*i) ,1 = 0 (1) 


Momentum 


3uj 1 

P + # jUj + — 


3t ' “YMl 2 


P, i - Pr 


f 1 

u i,jj + - U J,U 
3 


(2) 
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3T 7 - 1 3P 

P — + pU t T j T i j - Wh = 0 

3t 7 3t 


(3) 


>ecies 


3Y 

P — + P* i*,i “ *,ii + w = 0 
3t 


(4) 


State 

P - pT (5) 

where the commas denote partial derivatives, the repeated indices 

» 

imply svunming, Pr is the Prandtl number, and Y represents the 
fuel mass fraction. Note that only one out of three species 
equations (fuel, oxidizer, and product) is taken into account 
from the simple chemical reaction model [21]. The following 
characteristic parameters are used to render the above equations 


dimensionless : 




P ~ P / P 0 • 

P = ?*/P 0 * 

T = T*/T 0 * 


Ui = UjVVo* , 

t = 

Xi = xjvr 

(6) 

M b = v o*/a 0 * / 

h = h*/C p *T 0 * 

W = W*a*/V 0 * 2 



in which J?* is the flame thickness given by a*/v 0 *, with a* being 
the thermal diffusivity, M b represents the Mach number; k*, the 
thermal conductivity; a 0 *, the speed of sound; h*, the combustion 
heat release; and w*, the reaction rate, whose dimensionless form 



9 


is 


w 


BzT 8 f 


\ 


P 

T 


n 

Y n exp[-E/T] 


(7) 


with z denoting the oxidizer-fuel ratio; n, the order of chemical 
reaction; E, the activation energy given by E ■ E*/RT 0 *; and B, 
the rate constant. The superscript * represents dimensional 
quantities and subscript zero gives the mean value at the flame 
edge. 

The solid propellant is assumed to be homogeneous, with no 
condensed phase chemical reaction. The dimensionless form of the 
heat transfer equation in the solid phase is 

t 

3T 

B — - + UjT 8/1 - fT, #ll = 0 (8a) 


Furthermore, assuming that the heat transfer is one-dimensional 
gives, 

8T S 3T S 3 2 T 9 

B + r { = 0 (8b) 

at ay ay 2 

where 

B - p 8 Vpo* , 
r = r*/r* with 

Here, r denotes the burning rate at the solid surface and 
subscript s represents the solid phase. The decomposition 
process of the solid propellant at the surface is assumed to 


S = k s *Cp*/Jc*C s * 

3 * _ *„ * / * 

r ■ P 0 v o /Pa 
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follow an Arrhenius law; thus, 


r = exu 


— F 

r 1 

1 1 


^8 

' T s 

1 00 
[ 1^ 

1 



( 9 ) 


in which E 8 is the dimensionless surface activation energy, E s = 
E 8 */RT 0 * , and T 8 is the mean temperature at the surface. The 
solid-gas interface boundary conditions are determined by the 
dimensionless mass and energy balances across the interface, such 
that 


r 3T ' 

l 

' 3T ' 

k 3y ) 

♦ " ^ 

i ay / 


+ rL 


( 10 ) 


with S = k*/k 8 * and L = (H + * - H_*)/C p *T 0 *. Here, L is the 
latent heat of vaporization of the propellant and H* denotes the 
enthalpy changes. The subscripts + and - represent the gas and 
solid side at the interface, respectively. 

When a small pressure disturbance occurs in the combustion 
chamber, every field variable will be disturbed from its steady- 
state value and can be expressed in the form, 


F ■ F (0> + cF l 1 5 + < 2 F* 2) + . . . 


( 11 ) 


where F « {/>, Uj, T, Y, P) and e represents the perturbation 
parameter. The superscripts in the parentheses indicate the 
perturbation order. Furthermore, assuming sinusoidal fluctuation 
of pressure with time renders the variables in a different form: 

X 1,2,... 


pt I ) _ £< I > e i I tit 


( 12 ) 
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It is important to recognize that the sources in the second 
order consist of inhomogeneous terms that describe the 
nonlinearities in terms of the product of two first order 
variables. Considering the physical quantity of f < 1) leads to 
the separation of each inhomogeneous term into a time-independent 
term and a term oscillating at the frequency of the second 
harmonic, i.e., 

F l <n F 2 tl> = - ($! 1 1 ‘ 1 ’ + * 1 ’fcj 1 1 *e l 2wt ) (13) 

2 


with the dagger representing the complex conjugate. 

Consequently, the dependent variable F <2> may be rewritten as ^ 

f < 2> . $ r (2> + ^ D (2) e i 2 «t (14) 

For a higher order, the same argument is applicable. 

Substituting Eqs. (11) -(14) into Eqs. (l)-(5) and rearranging 
separately in the order of perturbation yield the final form of 
the governing equations corresponding to each order. 

Steadv-State Governing Equations 

The one-dimensional steady-state governing equations are 
given as follows: 

Continuity 

P ( 0 V°> = 1 (15) 


Energy 


3T 


( o ) 


2 m ( 0 ) 


a‘T 


9y‘ 


ay 


= w l0> h 


(16) 
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Species 

3Y 1 0 > 8 2 Y ,0> 

= -w* 0) (17) 

9y 9y 2 

State 

^ t o » T t o » _ 1 (18) 

Note that the uniform pressure is retained throughout the flame 
thickness (P t0> = 1), and from Eqs. (6) and (11)# both the 
dependent variables, p l0> and T ,0) , are equal to unity at the 
flame edge, which is far from the origin on the scale of H * . 

» 

From Eqs. (15) and (18), we have 

v ,0> = T ,0> (19) 

and from Eqs. (16) and (17), y <0> can be expressed as 
1 

Y (0> = - (1 - T (0> ) (20) 

h 


Therefore, 


w ,0> 



2 

exp [ -E/T ‘ 0 1 ] 


( 21 ) 


where 5=0 and the second order chemical reaction is assumed. 

Now, the equations are expressed in terms of temperature in 
the steady-state; thus, only the solution of Eq. (16) is 
required. Flandro [23] suggests a simple analytical model of the 
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mean flow field to facilitate multi-dimensional analysis in the 
higher order system. This model is given in the form 

u, t0> = u c [l - exp(-y/R c ) ] i + v ,0) ] (22) 


Here, u c describes the flow speed along the local streamline and 
R c is a dimensionless distance from the solid surface. The 
following boundary conditions are used in the steady-state 
solution: 

At the flame edge. 


Y ,0> ■ 0 (23) 

T ,0> = l. (24a) 

T 

or 

dT t0> 

= 0 (24b) 

dy 

At the solid phase, Eg. (8) gives 

T 8 t0> = (T s - T c ) e y/ * + T c (25) 


where T c is the propellant cold side temperature. The continuous 
temperature condition at the interface requires that 

T io> = T <o> = f (26) 

The matching condition across the interface can be obtained by 
substituting Eg. (25) into Eg. (10), which yields 


dT 10 ’ T 8 - T c 

= + L 

dy + n 


( 27 ) 
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In solving Eg. (16), note that a correct rate constant B in Eq. 
(21) has to be determined such that the system satisfies the 
boundary conditions at the flame edge as well as at the 
interface. 

Higher Order Governing Eguations 

The higher order governing eguations can be expressed in a 
single form because only the source terms are different from each 
other. Presenting the source terms on the right-hand side of the 
eguations in terms of G s , the governing eguations are represented 
as follows: 

Continuity 


ilwp 1 1 * + ( P <0 »Ui 


( i > 


* < i >„ 

P Uj 


( 0 > 


(28) 


Momentum 


4 Tt ,« { 0), A . <i> j. r« , o>, 1 <0> , A 1 ci) + ( 0 ) . Mi)* 0 ’ <o) n 

Hup Uj + [ p U^jUj + p U^jUj + p U^jUj ] 


’. ( I ) 


( 0 ) 


7Ml 


A ( I ) 

P,i - PH 


A ( I ) X A ( I ) 

u i,n + 7 

3 


- G 


2 i 


(29) 


Energy 


c o > 


ilwp* ® ’T* 1 * + [p (0> Uj 1 0 j + A p«‘-u 1 '«'T, i ] 


7-1 


- iiw 


A ( I ) 


p (I > _ _ £ci> h 


(30) 


Species 


iiw p * 0 * Y 1 1 ’ + [p l0, u i ,0> Y*j > + P 


< o ) , 


< I ) < 


( 0 ) 

/ i 


A ( I ) 
f> { 1 Uj 


( 0 ) 


( 0 ) 
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A I I ) A f » i 

- Y /U + V ,M = G t 


( 31 ) 


State 


£( i ) 


- P , 0 , T (I> - p tI, T (0> = G t 


( 32 ) 


and the reaction rate is given by 


w ,n = w ,0) 


( o > 


( i > 


E 2 

£( I ) + y< I ) 


,( 0)2 


( 0 > 


+ Gi 


( 33 ) 


Here, G 9 is given as follows: for the first order system (1= 1) , 
G = 0 ; for the second order system (I = 2 ) , 


r - _ * 1 > n < 1 > \ 

G 1 (p ) , l 

r - - / ^ A < 1 ) n 1 1 > 4 . « ( o ) A 1 1 * A ( 1 ) . A ( 1 ) „ 1 0 1 A ( 1 ) 

g 2 i “ -( 1 «P Uj + p U l; jUj + p U iy jUj 

+ P u i , j u j > 


g 3 = -{i«p ,l) T ll> + p (0, u 1 ,1, t;;' + p« i, u 1 ,o, t;i 


A ( 1 ) A 


( 1 )„ ( 0 ) 


A ( 1 ) 


. A ( 1 > A ( 1 ) m ‘ ° * s 

+ P u i T , i > 


A * • ^ A 


A ( 1 ) A 


A ( 1 ) 


g 4 = -{iwp ,i, Y ,i) + p ,o > + p ,i, u 1 ,o, Y; i 

( 0 ) 


4. A ( 1 > , A , ( l ) V \ 

+ P U i Y ,i > 


( 1 >m < 1 ) 


f 


w 


( 0 ) 


( 0 > 


2 E 


( 0 >m( 0 ) 2 


( 1 >m( X ) 


p‘ °'T 


p‘ 1 'T 


< 0 > v ( 0 > 


p Y 


A ( 1 ) v ( 1 ) 

p * 


E 

r E 1 


'Ip ( 1 ) 1 

2 2E 

. A ( 1 ) m ( 1 ) . 

ry‘ 1 >\ 

ip ( o ) 

l2T ,0> J 


T ( 0 ) 

v I } 

Y< o ) T ( 0 ) 2 

[ Y (0)J 


( 34 ) 
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Note that, including the pressure coupling, the velocity coupling 
is significant in the source terms. 

At the solid phase, Eg. (8) gives 

A aT 8 lI> * 3T 8 (0 » 3 2 t 8 1 1 ’ 

iluBTg ‘ 1 * + r (0) — + r' 1 ’ — £ = G 7 (35) 

dy dy dy 2 


with 


3y 

Linearizing Eg. (9) results in 

r lI> = r ,0, c l (T 8 ‘ I) + Gg ) (36) • 

where 



2 T 8 1 


Substituting Eg. (36) into Eg. (35) and solving analytically 
yield 


3 T‘ i } 

4- CL 

* I 

C lC 2 1 

-i i 

c t c 2 

i — T 

ay 

+ g 9 
+ / 

. % 

x ilwS y 

• C j 1 j 

ilwBC ^ 


in which 


C 3 C 2 

r i ] 

i 

C 4 ) 

C, = 

— X 2 + ““ 

+ - 

X 2 c 4 + C 5 - 

ilwSJ 


% 

k £ i 


and 
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1 

\ 1 = — [1 + (1 + i4IuBf) 1/2 ] 
2f 

Ci - e 8 /t 8 2 
1 _ 

c 2 - - (T, - T e ) 
f 


c, = c, 


Cl 

1 

\ 2 


$< x ) 2 


c 2 r 


( 1)2 


2u 2 B 2 f 


MX, “ X.) 


({X t 2 - Xj - i2 wfi 


( 1 ) 


T 8 


( l ) 


c,r 


( i ) 


iufi 


» 


Equation (37) gives a boundary condition at the surface; other 
conditions are as follows: 

At the flame edge, 


A 

Y 


( I > 


0 


(38) 


Assumption of an isentropic flow near the flame edge gives the 
temperature conditions. These conditions are equivalent to Eq. 
(3) after deleting the thermal diffusion and reaction terms. 
Noting that the steady-state temperature gradient at the flame 
edge is almost zero, this condition can be expressed in the form 


ilu p 


( 0 > m ( I ) 


„< 0 >„ 

p Uj 


( 0 >£ 


( I ) 
/ i 


- ilu 


7-1 


P lI > = G 


1 0 


(39) 


with 
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'10 


-{ iwp* 1 * T 


( l ) 


p ,0> u 


( l Jm* 1 * 
i T Z i 


; ci »ui ,o »T 


( i ) 

,i > 


Since the flux of each reactant species is always a fixed 
fraction of the total flux, the fuel mass flux fraction m f can be 
derived from Eq. (4) as 


pUiY,i -Y #li - 0 

For a one-dimensional expression, m f is given as 


(40) 


1 dY 

- Y (41) 

m dy 


where m is the mass flux equal to pv, and assumption of the 
constant burning rate at an instant has been used. At the 
surface, 


m f 


( i > 


A 1 dY (I> 1 A dY <0 > 

Y (1> + m lI> + G 

m <0> dy m <0)2 dy 


(42) 


in which 


1 

' m (1) 

dY * 1 * 

dY‘°> 

[ m ( 1 * | 

2 - 

m ,0 > 

m <0) 

dy 

dy 

m (0 > J 



Using the relationship m‘ 1 * = p s r ,I> yields 


A dY (I> dY ,0 > A 

Y ,n Ci (T 8 <n + G 12 ) (43) 

dy dy 


with 
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The normal velocity component is derived from the mass 
balance condition at the interface such that 



with 

1 3u ll> 1 a 2 u <0> 

G £ii> 

14 ilu& 3y + 2I 2 w 2 B 2 3y 2 

Note that G 7 - G l4 are valid only for 1=2. For higher order 

systems, Eqs. (37) -(39) and Eqs. (43) -(45) are used as boundary 

conditions to solve Eqs. (28)-(32). It is necessary to have more 

conditions for the density at both sides for better solutions, 

and the pressure fluctuation has to be forced at the flame edge. 

In the case of the second order time- independent system, care 

must be taken to use boundary equations (35) and (39) . 


£< i > 2 

+ 
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3. NUMERICAL METHOD 

Each set of governing equations, subject to appropriate 
boundary conditions, is solved using the Galerkin finite element 
method. The boundary equations are imbedded in the total matrix 
equations by means of Lagrange multiplier X Y , such that the 
overall global matrix has the form 

iI«A afi + B 

in which Xg represents the solution vector, X fl ** [p fi , u fil , Tg, 

Yg, Pg], and G a denotes the inhomogeneous source terms valid for 
the second order. The second row, q 7g Xg = b Y , represents the 
boundary equations, where 7 = l,2,...,m, m being the number of 
equations. The solution of Eq. (46) at a given frequency is 
obtained by imposing the Dirichlet condition of the pressure at 
the flame edge. Note that the first row of Eq. (46) gives the 
eigenvalue problem when G a = 0 and from which the natural 
frequency of the system is obtained. The finite element 
formulation contains two different kinds of test function used to 
represent the volume and surface integrals in the domain. The 
total number of field variables would be reduced by one if the 
density or pressure were eliminated using the perfect gas law. 
However, the stability of the matrix is doubtful. Before 
calculations, the following is expected: since the Mach number 

is generally very small, the coefficient of the pressure term in 
Eq. (2) dominates the system unless the frequency considered is 
large enough such that other coefficients containing the 
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frequency factor become comparable in magnitude. Therefore, in 
lower frequencies the pressure gradient has to be negligible, 
thus resulting in constant pressure. On the contrary, the 
gradient will become significant in higher frequencies; this 
results in pressure variance. In the latter case, severe 
pressure changes will occur if unbounded at the solid surface. 
Note that care must be exercised in expanding Eq. (9) due to the 
appearance of the exponential growth effect. For the steady- 
state case, as mentioned earlier, we calculate the eigenvalue B 
in which necessary initial conditions are satisfied. However, as 
will be discussed later, the eigenvalue is referred from the 
result of reference [21] for this study. 


4. DISCUSSION 

All the perturbed governing equations in the higher order 
systems having variable coefficients basically depend on the 
steady-state temperature distribution in the domain, as shown in 
Eqs. (15) -(18). Figure 1 shows a typical steady-state 
distribution of the field variables, including the reaction rate 
for an adiabatic flame with a second order chemical reaction 
mentioned in T'ien [21]. For verification, the following 
parameters are utilized: z = 1, 8 ■ 0, T c - 0.15, T s = 0.35, f = 

£ = 1, E = 10, E s = 4, L = 0.15, and B = 1000. The other 
parameters are taken from Flandro [23] and are as follows: 7 = 

1.2, M b = 0.003, u c = 1.0, R c = 5.0, Pr = 1.0, and h = 1.3. 
Representative dimensional parameters corresponding to the 
dimensionless values are given in Table 1. 
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Table 1 Typical value and range of parameters 


Parameter 

Typical 

Value 

Range 

Physical 

Variable 

Typical 

Value 

Pa 

1000 

250-1000 

-g/cm 3 

1.5 

p g 

1 

- 

g/cm 3 

1.5 x 10“ 

T c 

0.15 

- 

°K 

300 

t 9 

0.35 

- 

°K 

700 

Tf 

1.0 

- 

°K 

2000 

E 

10 

4-15 

cal/gmole 

40 X 10 3 

E s 

4 

2-10 

cal/gmole 

16 X 10 3 

C P 

- 

- 

cal/g°k 

0.33 

c 8 

- 

- 

cal/g°k 

0.33 

k g 

- 

- 

cal/cm°ksec 

1 

o 

H 

X 

in 

*s 

- 

- 

cal/cm°ksec 

5 X 10" 4 

m 0 

1 

- 

g/cm 2 sec 

0.4 

P 0 

1 

- 

atm 

9.5 

* 

2 

0.5-2 

- 

- 

Yf 

0.5 

0 
• 

1 

o 

• 

00 

U1 

- 

- 

U 

- 

10* 3 -10 2 

H 2 

- 

7 

1.4 

- 

- 


AH 

0.15 

0.05-0.3 

cal/gmole 

2.8 X 10 3 
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As mentioned earlier, the natual frequency of the given 
system is obtained from the homogeneous solution of the total 
matrix equation (46) without the boundary conditions. As 
previously suggested [28], the result shows that the active 
energy transfer between the acoustic wave and the combustion 
process occurs mostly in a lower frequency range which leads to 
acoustic instability. In this study, most of the frequencies are 
clustered in w < 20; hence, the frequency range of interest is 
chosen between u = 10' 3 and w » 10 2 and it extends up to u = 500. 

The thickness of the burning zone is assumed to be 
negligibly small compared with the wavelength of the acoustic 
oscillation; thus, the pressure is approximately uniform f 

A 

throughout the domain of study and varies only with time. From 
this point of view, one of the most significant aspects of the 
present study is the fact that the oscillating pressure is no 
longer taken to be uniform at any instant, but is regarded as a 
spatially nonhomogeneous time-dependent source term. 

Consequently, it allows us to investigate the response of a 
specific propellant at significantly high frequencies and to find 
the response in the long flame of a double-base propellant. The 
frequency limit has usually been determined by the reciprocal of 
the characteristic time in Eq. (16) . However, is is inversely 
proportional to the square of the burning rate; therefore, the 
limit cannot be constant, but varies with the pressure 
fluctuations. This argument is verified in the present study 
using two different cases; (l) increasing the order of 
perturbation decreases the upper limit of the frequency and (2) 
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increasing the flame thickness also decreases the upper limit. 
Flandro [23] defines the non-uniform pressure coefficient 
explicitly in terms of the position and incident angle outside 
the combustion region; however, further discussion concerning 
this subject is not available. 

First, we attempt to verify the results of the new model in 
the one-dimensional problem. The calculation is actually 
conducted in multiple dimensions, but the boundary conditions are 
chosen approximately as if the gas flow seems to act one- 
dimensionally. Data at the center nodes of the domain are used 
to evaluate the results. Figure 2 demonstrates distributions of 
component fluctuations in the first order at w ■ 1. The results 
are comparable to those of T'ien [21]. In Fig. 2, it is also 
shown that the highest temperature fluctuations arise at the 
point where the maximum reaction occurs, other field variables 
have their maximum/minimum values at the point where the gradient 
of the reaction is maximized. The temperature amplification at 
the surface changes the burning rate in Eg. (36), while the 
velocity at the flame edge represents the acoustic admittance, 
whose magnitude and sign indicate the instability of the system. 
Note that the pressure remains constant, implying that the 
acoustic wavelength is larger than the flame thickness in this 
case; consequently, the imaginary part of the velocity approaches 
a constant slope at the edge. 

Distributions of the field variables against frequency for 
the first order are calculated and shown in Figs. 3-7. The 
amplitude of the pressure disturbance playing the role of the 
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forcing function is set to unity at the - flame edge. The results 
show that, at certain lower frequencies (« < 10), the pressure is 
constant but varies in the higher region (w > 10) . Thus, the 
limit of the constant pressure assumption is clearly recovered. 
General trends show that the distribution profiles are divided 
roughly into two groups at about « = 1; one is for w < 1 and the 
other for « > 1. In the lower frequencies (quasi-steady region), 
the variables, keeping similar distribution profiles, change 
their magnitude negligibly along the frequency. This trend is 
also true for the higher frequencies up to w - 100. The only 
exception is at w * 1, and it gives very different distribution 
profiles among others. The figures also reveal that the results 
are closely related to the chemical reaction distribution. 

The overall density distributions based on a second order 
chemical reaction are depicted in Fig. 3. The changes in 
magnitude along the flame are more significant in higher 
frequencies than in lower frequencies. These changes seem to be 
directly related to the fuel species (Fig. 7) and implicitly 
related to the temperature (Fig. 6) . The imaginary parts of the 
density represent the phase shift from the incident wave, and 
these are almost zero except for u = 1. 

At the surface, positive magnitude implies a stagnant 
phenomena caused by decreasing the velocity, although the mass 
flux increases at a higher pressure level. Special attention is 
invited to the profile at « = 1, where the profile is entirely 
different from others and some portions of the flame have 
negative amplitudes. The mass balance predicts a faster velocity 



26 


at the negative portion, as shown in Fig. 4. It can be said that 
the combustion system is most sensitive to the pressure 
fluctuation at u - 1 in the frequency region considered. If the 
upper limit of the frequency is extended, the profile is 
reversed, but with a similar trend of periodicity. A simple 
chemical reaction model restricts a realistic discussion in 
detail since, for most propellants, it is more complex than it 
implies. 

The normal velocity distributions are shown in Fig. 4. Two 
kinds of profile are obvious: one with positive slope and one 
with negative slope. The latter contains most of the 
distributions in the lower frequency region, with some f 

exceptions. Note that a different profile appears at « = 1. 
Rearranging the real part of the velocity at the flame edge gives 
the distribution of the acoustic admittance, whose magnitude and 
sign indicate the amplifications or damping ability of the flame 
subject to the acoustic disturbance (Fig. 5) . Figure 5 reveals a 
resonance in the condensed phase near u = 0.01, indicating that 
the system is unstable. This verifies the early result of 
Denison and Baum [4]. Some negative peaks exist at the other 
frequencies, indicating that the resonance in the gas phase tends 
to damp the oscillatory motion. Figure 5 also shows the system 
to be unstable at most higher frequency regions. The real part 
of the burning rate in Eq. (36) gives a similar trend to the 
acoustic admittance at the quasi-steady region, although the 
magnitude is significantly different. But these trends differ 
from each other at the higher region, as indicated in [21] . 
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Thus, the burning rate is not representative as a stability 
measure for a higher oscillatory case. Over u = 100, the profile 
tends to have a second mode oscillation as the pressure varies in 
the flame zone. 

Figures 6 and 7 illustrate the temperature and fuel species. 
At the lower frequency region, changes of magnitude are 
insignificant while at the higher region, such changes become 
significant. Because the temperature increases with fuel 
consumption, the distribution curves are in opposing directions. 
It is also seen that linearly diminishing the fuel affects the 
temperature changes slowly. Note that the difference in the fuel 
amount at the surface implies the change in the burning rate f 
affected by the disturbances. Significant changes of variables 
are also given at w - 1. 

In the first order system, the constant pressure field is 
valid until w = 10; above that frequency the pressure varies. 

This result gives the limit of the uniform pressure assumption. 
Furthermore, it shows that up to u = 100, the magnitude grows 
linearly starting from the flame edge where the Dirichlet 
condition is imposed, and then begins to oscillate. 

As previously indicated, each variable of the second order 
response to acoustic disturbance has two components: one time- 
dependent component that oscillates at twice the fundamental 
frequency and one that is time- independent. The latter 
represents a shift in the mean value, thus causing a shift of the 
mean burning rate. The right-hand side of the second order total 
matrix equation consists of corresponding nonlinearities in terms 
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of the product of two first order variables. These 
nonlinearities function as source terms. The computations are 
performed for the time-independent component using the Dirichlet 
condition of the pressure to be unity at the flame edge. 
Distributions of the field variables as well as the shift of the 
burning rate are investigated. Figure 8 shows typical 
distribution profiles in the second order at « - 1? at this 
frequency, the constant pressure is retained. 

The variables follow trends similar to those of the first 
order, although the amplifications affected by the existence of 
nonlinearities in the higher order are different. The trends 
still show a small discrepancy at the flame edge as in the first 
order. Figures 9-12 illustrate the behavior of each variable 
against the frequency. The general tendency of the second order 
is to affect the flame toward stability in the lower frequency 
region. Note that the upper limit of the frequency for constant 
pressure assumption has to be reduced by one half. The 
computational results show that the pressure changes from u = 5, 
which is half of the limit frequency in the first order. Thus, 
the limit should be determined by considering the order of 
perturbation involved in the calculations. 

At u < 100, the variables have the same order of magnitude 
as rhat of the pressure imposed. However, the velocity changes 
significantly from w = 1 due to the pressure change from that 
frequency. Thus, a higher order effect may not be negligible 
unless the perturbation parameter e has an order of magnitude 
less than the reciprocal of the highest order in the problem. At 
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lower frequencies, as shown in the first order, the distribution 
profiles are very similar; at higher frequencies, they differ 
from each other. 

Finally, the burning rate offset is calculated and given in 
Fig. 13. The offset is relatively small in the quasi-steady 
region, but increases with oscillatory motion along the 
frequencies. It has a negative sign in the entire frequency 
range, indicating a decrease of the burning rate subject to the 
acosutic pressure oscillations. This verifies the experimental 
result of Price [39], and the averaged curve looks similar to the 
analytical results of Friedly and Petersen [10]. 

Parameter studies are conducted for the first order and t 
summarized as follows. Decreasing the density ratio fl affects 
the variables shifted slightly to the negative direction, keeping 
the distribution profiles constant. Changing the latent heat of 
solid L exerts a negligible effect on the variables, but a very 
small value of L shifts the system toward instability. 

Increasing the surface activation energy E $ or the gas phase 
activation energy E reduces the magnitude of the variables, 
keeping the same profiles. Changes of the rate constant and 
viscosity effect coefficient strongly affect the system, such 
that every aspect discussed herein will change. 

The present study could be extended to the multi-dimensional 
case by introducing the appropriate axial mean flow field [46]. 

It is well recognized that fluctuation of the gas velocity 
parallel to the propellant surface affects the burning rate in 
terms of velocity coupling; therefore, this quantity must be 
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considered together with pressure coupling for a satisfactory 
measure of stability. A simple calculation has been accomplished 
using the artificial axial flow velocity in Eq. (22) . However, 
difficulties of the boundary conditions could not be eliminated. 

A test run shows that the existence of a small amount of the 
axial flow reduces the range of dispersion of the variable 
distribution profile from each other in the frequency region that 
leads the system toward stability. 

5. CONCLUSION 

A multi-dimensional numerical model for the premixed flame 
acoustic instability is proposed and solved using the finite ^ 
element method. The governing equations are perturbed to the 
second order and formulated with Galerkin finite elements. The 
gaseous flame is assumed to be simple and homogeneous, and the 
Arrhenius manner of decomposition is implemented with no 
condensed phase reaction. The results have direct bearing on the 
validity of published theories of solid propellant combustion 
instability at the lower frequency region where the uniform 
pressure is valid. Extended studies are made on the higher 
frequency region and the results are discussed. Under the 
restricted boundary conditions, the following conclusions, based 
on numerical calculations, are reached: 

(1) The pressure is assumed to vary in the domain of study, and 
calculations based on nonuniform pressure indicate that for 
« > 10, there is a significant deviation from the uniform 
pressure assumption for the first order. 
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(2) For the second order system, such deviation occurs at a 
lower frequency which is half of the first order frequency 
limit. 

(3) Investigation of the distribution of variables shows that 
the acoustic instability is likely to be most critical at 
u « 1, while the acoustic admittance at the flame edge 
indicates a negative sign. 

(4) The oscillatory amplification or damping ability of the 
flame is recovered in the quasi-steady region and, at a 
higher frequency, moderate amplification effects are 
obtained. 

(5) The burning rate is directly related to the acoustic 
admittance only at the lower frequency region and its 
negative offset phenomena have been valid in the second 
order perturbation study. 

(6) Second order effects may cause the instability to be more 
critical in some cases and negligible in others. 

(7) Multi-dimensional instability calculations can be achieved 
using this model so far as a realistic mean flow field is 
clarified with proper boundary conditions. 
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Fig. 1 Steady-state distributions of field variables 

in flame zone. Parameters used in calculations 
are given in Table 1, with Pr = 1 and = 0.003 
Mean values at the flame edge are used to non- 
dimens ionalize the variables. The higher order 
calculations are based on this result. 
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Fig. 5 Acoustic admittance and burning rate vs. frequency. 
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Fig. 7 First order species ("fuel) distributions vs. frequency. 
















48 



Figv; 13 Offset in burning rate for the second order time- 
independent system. Calculations are from Eq. (36) 
using Figs. 6 and 11. The negative sign in all 
the frequency region indicates that the burning 
rate always decreases by the acoustic pressure 
oscillations reported by Price C39J. The average 
trend is comparable to the analytical result of 
Friedly and Petersen C103. 




